/aj -&> v - 
=^ 7 75 9*7 

■718. 

ADAPTIVE GRID EMBEDDING FOR THE TWO-DIMENSIONAL 
FLUX-SPLIT EULER EQUATIONS 


by 

Gary Patrick Warren 


A Thesis 

Submitted to the Faculty of 
Mississippi State University 
in Partial Fulfillment of the Requirements 
for the Degree of Master of Science 
in the Department of Aerospace Engineering 


Mississippi State, Mississippi 


May 1990 


(NASA-CR-186533) AHA^TTVE GRID EMBEDDING 
FOR THE TWO-DIMENSIONAL FLUX-SPl IT EULtR 
EQUATIONS M.S. Thesis (Mississippi State 
Univ.) 77 0 cscl - 12A 

G3/64 


N Q 0-21 571 


Unc 1 as 
0277537 



ADAPTIVE GRID EMBEDDING FOR THE TWO-DIMENSIONAL 


FLUX-SPLIT EULER EQUATIONS 


by 


Gary Patrick Warren 


APPROVED: 


Professor of Aerospace Director of Graduate Instruction 

Engineering (Major Professor) College of Engineering 


Associate Professor, Graduate 
Coordinator, Department of 
Aerospace Engineering 


Dean of the College of 
Engineering 
Aerospace Engineering 


Professor and Head of the 
Department of Aerospace 
Engineering 


Vice President for Graduate 
Studies and Research 


May, 1990 



ACKNOWLEDGMENTS 


The author wishes to thank the NASA Langley Research Center for their 
financial support for this project. Also, thanks go to Jim Thomas for helpful 
discussions and Jerry South for his patience. Special thanks to W. Kyle Anderson 
for always having the time to be a teacher and a friend. 


n 



ABSTRACT 


Gary Patrick Warren, Master of Science, 1990 

Major: Engineering, Department of Aerospace Engineering 

Title of Thesis: Adaptive Grid Embedding for the Two-Dimensional Flux-Split 

Euler Equations 

Directed by: Joe F. Thompson 

Pages in Thesis: 65 Words in Abstract: 185 

ABSTRACT 

l 2 

A numerical algorithm is presented for solving the two-dimensional flux-split 
Euler equations using a multigrid method with adaptive grid embedding. The 
method uses an unstructured data set along with a system of pointers for com- 
munication on the irregularly shaped grid topologies. An explicit two-stage time 
advancement scheme is implemented. A multigrid algorithm is used to provide grid 
level communication and to accelerate the convergence of the solution to steady 
state. Results are presented for a subcritical airfoil and a transonic airfoil with 
3 levels of adaptation. Comparisons are made with a structured upwind Euler 
code which uses the same flux integration techniques of the present algorithm. 
Good agreement is obtained with converged surface pressure coefficients. The lift 
coefficients of the adaptive code are within 2 \ a ( of the structured code for the sub- 
critical case and within ±\% of the structured code for the transonic case using 
approximately one-third the number of grid points. 
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NOMENCLATURE 


A,B flux Jacobians, and 

a speed of sound 

CFL Courant-Friedrichs-Lewy number 

C p pressure coefficient 

c chord length 

Cd drag coefficient 

ci lift coefficient 

e total energy 

F fluxes of mass, momentum, and energy 

Gi grid level i 

I ) -1 collection operator used for transferring information on grid level * 
to grid level t — 1 

IJ_ x interpolation operator used for transferring information on grid level 

i — 1 to grid level i 

A • < 

I t collection operator for residual 

A A 

ij unit vectors in x and y directions 

k n denotes ^ or r) 

i length of a face 

M Mach number 

Mach number in £ direction 

vi 


n 

unit normal 

h x ,hy 

x and y components of a unit normal 

p 

pressure 

Q 

conserved state vector, Q = \p,pu,pv,e\ T 

q« 

most current approximation to Q on grid level i 

7 

velocity magnitude 

R 

residual vector for mass, momentum, and energy 

S 

surface area 

8 

entropy, also parameter for flux limiter 

t 

time 

u,v 

contravariant velocities 

U,V 

cartesian velocities in x and y directions 

v, 

correction at grid level t 

x,y 

cartesian coordinates 

At 

time step 

a 

angle-of-attack 

P 

refinement parameter 

1 

ratio of specific heats, taken as 1.4 

V 

gradient operator 


accuracy parameter 

V 

undivided gradient operator 


general curvilinear coordinates 
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k parameter controlling spatial difference-scheme type 

p density 

t relative truncation error 

Subscripts: 

c information at coarse grid face 

/l, /2 information at fine grid faces 

oo denotes conditions at infinity 

Superscripts: 

n denotes time level 

+, — denotes direction from which information was obtained 


viii 


List of Figures 


Figure Page 

1 Arbitrary Control Volume 5 

2 Structured Grid 27 

3 Unstructured Grid 27 

4 Semi-unstructured Grid 27 

5 Single Grid Level Structure 28 

6 Multiple Grid Level Structure 28 

7a Cell-to-Face Pointer 29 

7b Cell-to-Node Pointer 29 

7c Face-to-Cell Pointer 29 

7d Face-to-Node Pointer 30 

7e Lower Link Pointer 30 

7f Upper Link Pointer 30 

8 Stability Characteristics for k = 1/3 31 

9 Stability Characteristics for k = -1 31 

10 Inflow/Ouflow Boundary 32 

1 1 Characteristics at Inflow/Outflow Boundary 32 

12 Embedded Ghost Cells 33 

13 Computational Stencil for Fine Grid Fluxes 34 

14 Computational Stencil for Coarse Grid Flux 34 

15 Multigrid V-cycle 35 

16 O-Grid Topology 41 

17 Voids in Grid 42 


IX 



18 Grid Remapping Criteria t 42 

19 Grid Adaptation Steps 43 

20 Initial Grid 47 

21 Initial Grid Pressure Distribution, = 0.63, a = 2.0° 47 

22 Initial Grid Mach Contours, = 0.63, a = 2.0°, AM = 0.05 48 

23 Leading Edge Initial Grid Mach Contours, Moo = 0.63, a = 2.0° 48 

24 2 Level Adapted Grid (V/?), M^ = 0.63, a = 2.0° 49 

25 2 Level (Vp) Pressure Distribution, M^ = 0.63, a = 2.0° 49 

26 3 Level Adapted Grid (Vp), Mqo = 0.63, a = 2.0° 50 

27 3 Level (Vp) Pressure Distribution, M^ = 0.63, a = 2.0° 50 

28 4 Level Adapted Grid (Vp), M^ = 0.63, a = 2.0° .51 

29 4 Level (Vp) Pressure Distribution, Mqo = 0.63, a = 2.0° 51 

30 2 Level Adapted Grid (Vq), = 0.63, a = 2.0° 52 

31 2 Level (Vq) Pressure Distribution, M^ = 0.63, a = 2.0° 52 

32 3 Level Adapted Grid (Vg), M^ = 0.63, a — 2.0° 53 

33 3 Level (Vq) Pressure Distribution, M^ = 0.63, a = 2.0° 53 

34 4 Level Adapted Grid (Vq), M^ = 0.63, a = 2.0° 54 

35 4 Level (V^) Pressure Distribution, M^ = 0.63, a = 2.0° 54 

36 4 Level (Vp) Mach Contours, M^ = 0.63, a = 2.0°, AM = 0.05 55 

37 Leading Edge 4 Level (Vp) Mach Contours, Mqo = 0.63, a = 2.0° .... 55 

38 Initial Grid Mach Contours, M,*, = 0.80, a = 1.25° 57 

39 Initial Grid Pressure Distribution, M^ = 0.80, a = 1.25° 57 

40 4 Level Adapted Grid (Vq), M^ = 0.80, a = 1.25° 58 

41 4 Level (Vq) Pressure Distribution, M TO = 0.80, a = 1.25° 58 

42 4 Level (Vq) Mach Contours, Mqo = 0.80, a = 1.25° 59 

43 Typical Quadrilateral 64 

44 Quadrilateral in a Uniform Cartesian Grid 64 


X 



List of Tables 


Table Page 

4.1 Expected effectiveness of various refinement parameters for inviscid, 

transonic flows over airfoil-like bodies 37 

5.1 Results Comparison for = 0.63, a = 2.0° 56 

5.2 Results Comparison for = 0.80, a = 1.25° 59 


XI 


CHAPTER I 


INTRODUCTION 

In spite of many recent developments in computational fluid dynamics al- 
gorithms and substantial increases in computer speed and memory, many flow 
situations still impose unsatisfactory requirements on computer time. These re- 
quirements are driven by the large number of grid points required to sufficiently 
resolve important flowfield features. Since most algorithms use structured grids 
without embedding, adequate resolution of important flowfield features requires 
either highly stretched grids or the inclusion of many unnecessary grid points away 
from the main areas of interest. 

As a solution to this problem, several adaptive techniques have been developed 
in recent years which attempt to concentrate grid points only in regions of high 
gradients. These techniques can be grouped into one of two categories. In the first 
category, higher resolution of specific areas is accomplished by redistributing the 
grid points of an existing structured grid into the regions of interest. The primary 
advantage of this method is its relative ease of implementation into existing flow 
solvers because no modifications to the data structure is required. Unfortunately, 
this technique often results in highly stretched and skewed grid lines and can leave 
some regions of the flowfield with inadequate resolution. Another method some- 
times used to enhance grid resolution is the embedding of grid points into specific 
areas of an existing grid based on some reasonable criteria such as high flow gra- 
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clients or high local errors. This allows improved resolution of important features 
with a relatively small increase in grid points and without adversely affecting 
other portions of the grid. Grid embedding techniques, however, are difficult to 
implement into existing computer codes because of the necessary data structure 
required for inter-grid transfer of information. The difficulties encountered in grid 
embedding, however, are more than compensated by increased capability. 

There have been several recent investigations into grid embedding using 
central-difference algorithms 1 ’ 2 ’ 3 . Berger 1 developed a grid embedding technique 
for a structured grid algorithm using a central-difference method developed by 
Jameson. This algorithm adapted to an estimate of the truncation error. More 
recently, Dannenlioffer 2,3 developed an adaptive grid embedding method that was 
used in conjunction with a central differenced node based scheme. Dannenhoffer’s 
method used an unstructured grid along with a multigrid algorithm. 

In recent years, upwind differencing has gained in popularity because of in- 
creased accuracy and robustness over its central differencing counterparts. The 
finite volume technique, as described in references 4 and 5, has been applied to a 
variety of flow problems. This scheme uses the flux- vector splitting of van Leer 6 or 
the flux- difference- split ting of Roe 7 coupled with an implicit, multigrid algorithm. 
This multigrid scheme is a well proven and documented algorithm which has been 
used mainly on block-structured grids with limited use on embedded grids 8 . 

In the present study, the same spatial differencing used in reference 4 is applied 
to a data structure designed for the use of an unstructured adaptive grid embed- 
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ding algorithm. Because of the unstructured grids, an implicit time advancing 
method is not easily implemented and so the equations are advanced in time with 
an explicit, two-stage time marching scheme. Since explicit time marching often 
results in poor convergence rates, a multigrid algorithm is employed. Use of the 
multigrid algorithm not only improves the convergence rate, but as importantly, 
it also provides a mechanism for inter-level grid communication. This thesis will 
present algorithmic details and results of this method. 
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CHAPTER II 


GOVERNING EQUATIONS 


2.1 Euler Equations 


The governing equations are the time dependent Euler equations, which ex- 
press the conservation of mass, momentum, and energy for an inviscid gas. The 
equations are given by 

^ + VF = 0 (2.1) 

ot 

where the state vector Q and the flux vectors F are given as 


r p 



L e J 


F = Fi i + F 2 j 



pu 


■ pv ' 

F, - 

pu 2 4 - p 

f 2 = 

puv 

puv 


pv + p 


- (e + p)u. 


-( e -f p)v. 


The equations are closed with the equation of state for a perfect gas 

P = (7 - 1)[? - /5(" ! + f ,J )/2] 


( 2 . 2 ) 


(2.3) 


(2.4) 


(2.5) 


The variables in equations 2. 1-2.4 can be non-dimensionalized by introduction of 
the following definitions: 


4 


V — 


V 


X 


p = 


p 


Poo 



* = v 


icioo 

T~ 


u 





a — 



( 2 . 6 ) 


The non-dimensional equations are not repeated here since introduction of 
the previous variables into equations (2.1)-(2.4) results in an identical form of 
equations as given above but with the ( ~ ) removed. 


2.2 Integral Formulation 

The Euler equations can be integrated over an arbitrary control volume, as 
shown in figure 1, with the resulting equations given by 

IIs^ dS+ IL {vf)ds=o (2 - 7) 



Figure 1. Arbitrary Control Volume 
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The first integral in equation (2.7) is then written as 


where Q is given by 



( 2 . 8 ) 


Q = SQ 


(2.9) 


and the average value of the state vector over the control volume is defined as 


*=Ws qds 


( 2 . 10 ) 


The second integral in equation (2.1) is converted to a line integral by use of 
the divergence theorem 


J j (V • ¥)dS = j>W Ad( 


( 2 . 11 ) 


where h is the outward pointing unit vector normal to the control volume. 


n = {.i + (ti 


( 2 . 12 ) 


In tliis form, the fluxes normal to the surface of the control volume are written as 


F = F • n 


pu . 

puu + i z p 

pl'V + lyP 

( e + p)U J 


(2.13) 


where U is the velocity in the direction of the outward pointing normal 


V — + £y v 


(2.14) 


6 



2.3 Flux- Vector Splitting 


The time dependent Euler equations form a hyperbolic system of equations 
and hence information is propagated along characteristics. Therefore, to correctly 
model the flow physics, some form of upwind differencing should be used in order to 
model the characteristic nature of the governing equations. The upwind method 
used in the current study is the flux- vector-splitting method in which the flux 
vectors are split into forward and backward running contributions and differenced 
accordingly. 

The split flux- vectors used in the current study were developed by van Leer 6 . 
These have been generalized to curvilinear coordinates by Anderson 4 and have 
since been widely used for flow computations on structured grids 4 ’ 5 ’ 8 . They are 
continuously differentiable at eigenvalue sign changes and allow shocks to be cap- 
tured with at most two interior zones. In practice, only one zone is usually 
observed 4 . After splitting, the flux vectors normal to the cell face given in equation 
(2.13) can be re-expressed as 

F ~ F + + F“ (2.15) 

where F + and F _ are given below for | |< 1 


F± 


/± 

/ ± 


771 Cl 3 3 

mass 



+ «] 

+ *>] 


(2.16) 
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where 


f ± m ... = ± l) ! /4 (2.17) 

=/ ± m..,({-(7 - l)tl 2 ± 2(7 - l)Pa 

(2.18) 

+ 2a 2 }/(-, J -l) + („ ! + „ 2 )/2] 

Here, ^ is the Mach number normal to the surface of the control volume. 

For supersonic flow (i.e.| M { |> 1), the fluxes are given as 


F+ 

= F 

F~ = 0 

(Mf. > +1) 

(2.19a) 

F" 

= F 

F + = 0 

(M<: < -1) 

(2.196) 


where F is given by equation (2.13). 
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CHAPTER III 


NUMERICAL METHODS 


3.1 Data Structures 

The governing equations are solved on a discrete curvilinear grid which con- 
sists of cells, faces, and nodes. Generally, a structured grid consists of a collection 
of cells which are organized into a single or multiple block structure as shown in 
figure 2. In each block an implied data structure exists in which the relation- 
ships between the nodes, cells, and faces Eire determined by simple addition or 
subtraction of integer values from the current cell index (i.e. i + 1, i — 1, etc.). 
Unstructured grid algorithms, however, rely on a system of pointers to determine 
node, cell, and face relationships needed to transfer information within the grid. 
Every cell, face, and node is assigned an integer value, as shown in figure 3, and 
access to each is obtained through the use of various pointers. For example, to 
obtain information stored at the nodes which define a cell, an array called a cell- 
to-node pointer must be accessed. The number and type of pointers needed is 
highly dependent on the computational stencil for a particular scheme. 

To efficiently employ adaptive grid embedding, it is desirable to remove the 
restrictions of the implied data structure required for single and multiple block 
grids. In the current work, the grids consist of a collection of cells which are aligned 
with the £ and ?/ coordinates but no block structure is required. In addition, 
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all information in the grid is obtained through a system of pointers as in an 
unstructured grid. The use of the £ and r/ coordinate system makes it simple 
to implement the spatial differencing techniques of structured algorithms in an 
unstructured framework. This approach can be referred to as a semi-unstructured 
method since some structure is required for the cell alignment but not the cell 
collection as shown in figure 4. 

There are two basic methods in which the embedded grid regions can be 
handled. The first method is to treat the entire grid as a single level and have 
special data structures to handle the hanging nodes as shown in figure 5. This 
can be somewhat tedious because some cells may have more than 4 faces. For 
example, cell a in figure 5 consists of 6 faces and 6 nodes. The second method, 
which is used in the current work, is to treat the embedded grids as separate 
levels as shown in figure 6. This allows every cell to have only 4 faces and nodes. 
This multi-level grid structure also facilitates the use of a multigrid algorithm to 
enhance the convergence rate of the flow solver. The use of upwind differencing, 
adaptive embedding, and multigrid, is accommodated with the following set of 
pointers: 

1. Cell-to-Face pointer - relates a cell to its four faces. 

2. Cell-to-Node pointer - relates a cell to its four nodes. 

3. Face-to-Cell pointer - relates a face to the cells used in determining quantities 
on the face. 

4. Face-to-Node pointer - relates a face to its two nodes. 
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5. Link pointer - relates a cell to the four cells below it and the four cells above 
it. The coarse cell which lies below four finer cells is referred to as a parent 
and the finer cells are referred to as children. 

These pointers, and the numbering system for each, are demonstrated in figures 
7a-f. For setting the boundary conditions, additional information is required which 
contains the face numbers of both the solid walls and the inflow/outflow bound- 
aries. The indexing used in the algorithm descriptions will be based on the use of 
these pointers. 

3.2 Flux Integration 

The equations are solved using a finite volume formulation which entails solv- 
ing for the average quantities in each cell. The flux integral in equation (2.11) is 
approximated by: 

/(F )di = (3-1) 

71=1 

where the subscript n is the face number (see figure 7a) and is the length of 
each face. Note that the flux balance given by equation (3.1) is based on outward 
pointing normals whereas the normals for the semi-unstructured grid have been 
chosen, for convenience, to always point in the positive £ and ?/ direction. The 
resulting flux integration (with splitting) for a cell is written as 

/(?)<*{ = - [Fj'(Qf) + Fr(Q+)K, + IF+(Qj ) + F 2 -(Q +)]( 2 

J( (3.2) 

+ [F+(Q J -) + F,-(Q+)K S - IFJ-(Q 4 -) + F 4 -(Q +)]( 4 
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The fluxes F* are constructed using the same equations as but normals point- 
ing in the positive coordinate directions are used in place of outward pointing 
normals. Here, £ n is the physical length of face n and Q* are the state variables 
interpolated to the faces from the positive and negative directions. The interpo- 
lation formulas, using the stencil shown in figure 7c, are given by 

Q/ ace (n) ~ QceH(2) + ^[(1 “* + (1 + «)A + ] QceU(2) (3.3) 

Q/ ace (n) = Qcell( 3) — £^[(1 + + (1 “ *)A + ]Q C eU(3) (3.4) 

where A + and A- are forward and backward difference operators defined as 

A+Q cell(n) = Qcei?(n+1) “ Q cell(n) (3.5) 

A-Q cell(n) = Q cell(n) ~ QceU(n-l) (3-6) 

For k = — 1 , these interpolation formulas give a fully one-sided extrapolation of 
the dependent variables to the cell faces which results in a second-order accurate 
difference formula for equation (3.2) on a uniformly spaced grid 5 . Similarly, the 
use of k = | results in a three-point interpolation formula which produces a third- 
order accurate difference approximation to equation (3.2) on a uniform grid. 

When using a three-point interpolation formula, oscillations in the solution 
can occur in regions of strong gradients such as shocks. To eliminate these oscil- 
lations, flux-limited interpolations 5 are employed which are given by 

Q facc(n ) = QccU(2) + “ «^)A_ + (1 + **) A + ]Q C eU(2) (3.7) 
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(3.8) 


Q/ a «(„) = Q«n(j) - j^l(i + «»)A_ + (i - ««)^+]Qc«h(j) 

where 

2 A A _ -}- c 

9= (A + )2 + (A_)2+ £ 
and e is a small number (e = 10~ 6 ) preventing division by zero in regions of null 
gradients. The variable <j> is a user defined parameter which is equal to 1 but can 
be set to 0 for a one-sided first-order method. 

3.3 Time Advancement 

The solution of equation (3.2) is advanced in time using the modified Euler 
method 9 , which is an explicit 2-stage scheme of the form 


Q* - Q n + A*R n (3.10) 

Q n+1 = ^(Q n + Q*) + ^A<R* (3.11) 

£ it 

where the residual is given as 

R = -{ - [#+(QD + f r(<tf >]<1 + FJW) + f 2 "(Q 2 + )]< 2 

(3.12) 

+ [F+(Q,-) + f r(Qj-)]/» - [F 4 + (Q 4 -) + f 4 -(q+)]< 4 } 

The superscripts for the residual designate the time level of the state variables 
used in the flux evaluation. 

The solution for each cell is advanced at its own local time step corresponding 
to a CFL number given by: 

CFL = A#{| U\ + \ V | +a[\ gradtf) | + | grad(rj) |]} (3.13) 
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The contravariant velocity in the £ direction, U, is given in equation (2.14). The 
contravariant velocity, V , is obtained simply by replacing £ with rj. 

3.4 Stability Analysis 


The stability analysis follows from the work of reference 4. The stability 
characteristics are examined using the linearized two-dimensional Euler equations 
given by 


dQ 

dt 


-P( Q) 


(3.14) 


where 


dQ dQ 

F(q ) = A ^ +B W 


(3.15) 


The flux Jacobians, A = and B = are computed using the split flux 

vectors for the x and y directions respectively. The differentials are computed 
across the cell using the interpolation formulas given by equations (3.3) - (3.6). 
The Fourier mode for Q is assumed to be 


Q = Qe i0x e iry 


(3.16) 


where Q is an initial state vector. Substituting — P(Q) for the residual in equations 
(3.10) and (3.11) and replacing Q with its Fourier mode, a system of equations is 
obtained of the form 

Q n+1 = GQ" (3.17) 

The magnitude of the maximum eigenvalue of the amplification matrix G must 
be less than unity for stability. A computer program was used to cycle through a 


14 



fixed number of each of the spatial frequencies in the range 


0 < f3Ax,TAy < 2n (3.18) 

for a series of CFL numbers between 0.1 and 1.2. The generalized eigenvalue 
problem is solved using an IMSL 10 routine. In addition, the high frequency error 
damping properties (smoothing) of the scheme can be studied by examining the 
maximum eigenvalue in the frequency range between j and . 

Shown in figures 8 and 9 are the stability characteristics for k = 1/3 and 
k = -1 for a Mach number of 0.8 and a = 0.0. Both schemes exhibit poor high 
frequency error damping properties. Because the k = 1/3 scheme allows a higher 
CFL number along with higher spatial accuracy, this was the scheme chosen for 
the computations presented later. 

3.5 Boundary Conditions 

For setting the boundary conditions, there are three possible boundary types 
which must be considered. These include inflow/outflow, solid wall, and interfaces 
between embedded grid regions. All of the boundary conditions are implemented 
using “ghost” cells which lie outside of the computational domain. The solid wall 
and inflow /out flow boundary faces have only one ghost cell and hence the general 
four point stencil shown in figure 7c does not exist. The differencing on these faces 
is accomplished by setting Q + or Q“ (depending on the boundary) equal to that 
of the ghost cell. For the solid wall and embedded interface boundary conditions, 
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the ghost cells are updated before each time step whereas, for the inflow/outflow 
boundaries, the ghost cells are set to freestream values. 


3.5.1 Inflow/Outflow Boundary Conditions 

The inflow /out flow boundary conditions Eire developed using characteristic 
theory 11 which is summarized here. A local coordinate system (sc, y) is constructed 
orthogonal to the boundary as shown in figure 10. The characteristic equations 
are constructed assuming the tangential derivatives are negligible so that a one- 
dimensional analysis normal to the boundary can be employed. The resulting 
characteristics are sketched in figure 11. Assuming locally isentropic flow, the 
characteristic equations are given by 


dt 


ds 

— =0 

dt 

along 

C° = u 

(3.19a) 

dv 

— =0 

dt 

along 

C° = u 

(3.196) 

& 

H- 

II 

o 

along 

C ± = u ± a 

(3.19c) 


where the Riemann invariants, are given by 


fl±sfl± F^j < 3 ' 

The variables at the new time level /" + 1 are obtained from the Riemann invariants. 
Assuming that x points away from the computational domain, i? + and R~ are 
evaluated using the state variables from outside and inside the computational 
domain respectively. The normal velocity and speed of sound on the inflow/outflow 
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boundary axe obtained by adding and subtracting the Riemann invariants. 


Uface =\(R + + R~) (3.21) 

a face = {j -^-(R + -R-) (3-22) 

Once the normal velocity and speed of sound are calculated, the Cartesian 
velocities are determined by decomposing the normal and tangential velocity vec- 
tors 


^/ace — ^re/ face ^ref) (3.23) 

^ f ace — *Ve/ ^y(^/ace ^re/) (3.24) 


where the subscript ref represents freestreain values for inflow or from the cell 
inside the domain adjacent to the boundary for outflow. It should be noted that 
since the local coordinate x points away from the computational domain, inflow 
corresponds to u < 0 and outflow corresponds to u > 0. The variables h x and n y 
are the x and y components of the outward pointing unit normal for the face. 

The entropy s is determined using the value from outside the domain for 
inflow and from inside the domain for outflow. Once the entropy on the boundary 
is known, the density on the face is calculated from the entropy and speed of sound 
on the face 


r 


P face — 


face 

m lf & face . 


(3.25) 


The energy on the face is then calculated using the equation of state. 
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For supersonic conditions, the state variables are extrapolated to the bound- 
ary from the exterior for inflow boundaries and extrapolated from the interior for 
outflow boundaries. 

3.5.2 Solid Wall Boundary Conditions 

The variables on solid wall boundaries are computed by extrapolating from the 
interior of the computational domain. For the cases presented here, the pressure 
and density on the face are simply set to the values of the cell above it. The 
normal velocity is then set to zero so that the resulting Cartesian velocities on the 
face are given by 


u face — w cell ^x(^ceJf) 

(3.26) 

Vface ~Vcell ^y(^ceif) 

(3.27) 


Here, u C eil is constructed with a normal facing in the positive £ or r] direction and 
h x and h y are the x and y components of this normal. The density, pressure, and 
Cartesian velocities are then extrapolated to the ghost cells using 

Aghast = 2 <f>face “ <!>cell (3.28) 

Note that only the Cartesian velocities need to be extrapolated since the pressure 
and density in the ghost cell are the same as in the computational cell above it. 

3.5.3 Embedded Interface Boundary Conditions 

A difficulty encQuntered when using embedded grids is in determining the 
fluxes on the interface between the coarse and fine grids inside the computational 
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domain so that the accuracy and conservation properties of the flow solver are 
maintained. The method for obtaining the flux on the coarse grid is the same as 
that for obtaining the other coarse grid faces since the four-point computational 
stencil is complete. Obtaining the flux on the two finer faces can be accomplished 
several ways. The method developed here is similar to that of Berger 1 . A set 
of ghost cells is constructed outside of the interface region as shown in figure 12. 
Note that with the addition of the ghost cells, the four-point stencil needed for 
the interpolations to obtain the state variables on the faces is preserved, hi this 
region, the state variables are determined by bilinear interpolation from the next 
coarsest grid level using the most current solution on that level. This produces 
the interpolation stencil shown in figure 13 for the fine grid fluxes since all nine of 
the shaded coarse grid values are used to obtain the values in the ghost cells. The 
coarse grid flux stencil is shown in figure 14. 

This method of calculating the fluxes, however, does not enforce conservation 
at the interface region since there is nothing to ensure the sum of the fluxes on 
the finer grid equals the flux on the coarse grid. To enforce conservation at the 
interface, the fluxes from the two finer faces are added and injected onto the coarse 
grid before the coarse grid update 

t c t c =tfi( f i +*f 2 e f2 (3.29) 


3.6 Multigrid Algorithm 

The algorithm presented thus far is still incomplete because no communication 
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exists between the grid levels. The increased accuracy of the fine grids is not 
passed to the coarser grids, so the solution on the coarse grid remains unchanged. 
This is unacceptable since the ghost cells for the finer grids are determined by 
interpolation from the coarse grid solution. A multigrid algorithm is a natural 
way to obtain grid-level communication. The truncation error between a fine and 
a coarse grid is passed to the coarse grid in such a way that the coarse grid is 
driven by the fine grid residual. Hence, the fine grid accuracy is maintained on 
the coarser grid. 

The multigrid algorithm used in the current study is the Full- Approximation 
Scheme (FAS) which has been primarily used in solutions of the Euler equations 
as a method of convergence acceleration such as described in reference 4. Because 
of the embedded grids used in the current algorithm, slight modifications to the 
algorithm presented in reference 4 must be made since not all cells on the coarser 
grids have finer cells above them. 

The application of multigrid algorithms to Euler solvers has primarily been 
due to the convergence acceleration which is attained. The improved conver- 
gence is attributed to the fact that the multigrid algorithm efficiently damps low- 
frequency errors on the finer grids by using a sequence of grid levels denoted by 
Gat, Gn- i? where the coarsest grid level corresponds to Gj. Since the 

lower levels have greater spacing between grid lines, the low frequency errors on 
the higher levels appear as high frequency errors on the lower levels where they 
can be efficiently damped. 
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Explanation of the multigrid process is facilitated by first examining the Euler 
equations written in operator notation and discretized on a given grid, G jy, which 
is defined to be the finest level. 


I'aKQn) = 0 (3.30) 

where Q/v is the exact solution to the discretized system and L/y is the discrete 
steady state operator given by 

tw(Qw) = -{ - (Fj'(Qr) + *7(Ql')Ki + [P 2 + (Qn + iT (<#)]<* 

(3.31) 

+ [Pj + (QD + f,-(Q, + )Kj - [F, + (Q.~) + f.-(Q< + )K<} 

When solving iteratively, equation (3.31) is solved approximately at each time 
step so that the right hand side of (3.30) is not identically zero. This can be 
represented by introduction of a residual defined implicitly by 

^ n {^ n ) = (3.32) 

where is the most current approximation to Q;v and R;v is the residual. Note 
that the residual will be zero only when the approximate solution (q^) is equal 
to the exact solution (Qn)* 

Since it is the errors that are damped on any given grid, it is desirable to cast 
the governing equations for a given grid in terms of the errors. This is achieved 
by subtraction of equation (3.32) from equation (3.31) 

L/v(Q/v) - hTv(q^) = — R;v (3.33) 
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If the time stepping method adequately eliminates the high-frequency errors on 
the current grid level Gn, the residual equation, (3.31), may be adequately ap- 
proximated on the coarser parent cells Gn-\ by 

Ln-i(Qn-i) = ^’(-Rn) + L N _,(lX _1 q^) (3.34) 

Here I^ _1 is a volume- weighted collection operator that transfers dependent vari- 
ables from the children to the parents so that conservation is maintained 4 and is 
given by 

1 N lin = (3.J5) 

Here, V is the volume of each child cell. Similarly lj^ _I is the collection operator 

for the residual which is defined as 

i^- 1 R N = J]R A r (3.36) 

For equations (3.35) and (3.36) the summations are over the four fine grid cells 
(children) which make up the coarser parent cell. In the current study, the depen- 
dent variables in the ghost cells for the finer grid are not restricted down to the 
coarser grid. 

Note that equation (3.34) can be written as 

L/v-i (Q/v-i ) = Tjv-i (3.37) 

where 

tn-1 — Liv-i(ljv _1 qAr) — ijv~ 1 (R' v ) (3-38) 
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is the relative truncation error between grid levels G/v and Gn-i- For coarse grid 
cells with no finer cells above them, r/v-i is set to zero. The inclusion of j 
allows the solution in the parent cells to be driven by the finer grid so that the 
order of accuracy on the fine grid is maintained. The new residual which can be 
written in a form applicable to all grid levels, Gi is given by 

Rt = Li(q^) — T{ (3.39) 

This residual is now used for updating the solution on the current grid using the 
time advancing algorithm, equations (3.10) and (3.11), described earlier. 

To proceed down to the next coarsest level, G/v- 2> the analogous equation to 
equation (3.37) is given by 


LiV-2(QiV-2) = Tjy — 2 (3.40) 

where 

t n-2 = L N _ 2 (I^:^q^_ 1 ) - (3.41) 

Note that R;\r_i contains the relative truncation error between levels N and N — 1 
so that the relative truncation error on this grid is the sum of the truncation errors 
between levels N and N — 1 as well as N — 1 and N — 2. 

Thus far, the relative triuication errors between the parents and children 
have been succesively passed down to each of the parent cells so that the order 
of accuracy of the children is preserved. In this manner, the accuracy of the finer 
levels “drives” the accuracy of the coarser levels. 
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Since the lower levels damp out the low frequency errors of the higher levels, 
a correction from the coarse grid must be passed up to the higher levels. This 
correction V is formed on grid level N — 2 as 


Vn— 2 = q e N— , - Ijfrfow— , (s-«) 

This correction is then added to the dependent variable in the children cells 

<&-, - q c „-, + (3.43) 

where I^v — 2 1S an interpolation operator for passing information from a lower level 
to a higher level as described in Appendix A. A correction is then formed on level 
N — X and passed up to level N using equations (3.42) and (3.43) by replacing 
N — 1 with N and N — 2 with N — 1. Note that no additional iterations are 
performed on grid N — 1 before interpolation of corrections to grid level N. 

The process described above is a three-level “V”-cycle algorithm and is de- 
picted in figure 15. For further clarity, detailed steps for a three-level V-cycle are 
given below. 

1. Start on the highest grid level and advance the solution in time using equation 
(3.39) and r/v — 0. 

2. Calculate the residual with the most current values 011 the highest level from 
equation (3.39) and r/v = 0. 

3. For cells on G/v-i with children, collect the dependent variables from G/v 
using equation (3.35). 
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4. For cells on G/v-i with children, collect the residual from Gat using equation 
(3.36) and calculate the relative truncation error between the grids using 
equation (3.38). For cells on Gn-i with no children, set r;v-i = 0. 

5. Advance the solution in time on level G^r-i using equation (3.39) to calcu- 
late the residual. Since r^y-i has been previously calculated, the residual is 
calculated by evaluating LAr-i(q5v-i) from the most current values of the 
dependent variables on tills level and subtracting Tat-i* Again, note that 
tat-i — 0 for cells on Gat-j with no children. 

6. Calculate the residual on this level using equation (3.39). 

Rat- 1 = JjA r -l(qA/'-i) — r N - 1 

7. For cells on Gat -2 with children, collect the dependent variables from Gat-i 
using equation (3.35). 

8. For cells on Gat -2 with children, collect the residual from G w-i using equation 
(3.36) and calculate the relative truncation error using equation (3.41). For 
cells on Gat -2 with no children, set tat -2 = 0. 

9. Advance the solution in time on this level using equation (3.39) to calculate 
the residual on Gat- 2 • Since this is the lowest level used in the present exam- 
ple, several solution updates are performed to get an approximation to Qat- 2- 
During each step, the residual is updated to use the most current values of 
the dependent variables in La^^a^)- Note that r^v - 2 will not change. 

9. Calculate the correction on this level to give 

jN — 2 c 

N-2 - qiV-2 - jV — l^iV— 1 
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10. Pass the correction to the next finest mesh using bilinear interpolation and 


update the solution to give 

Qn + I/v-iVjv-i 

12. Calculate the correction on the N — 1 level with 

V»-1 = q e w _, - ijJ-’qS, 

13. Pass this correction to the highest level and update the solution 

<1 N + I^-iVn-1 

14. Advance the solution in time to smooth the errors. 
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Figure 2. Structured Grid 



Figure 3. Unstructured Grid 



Figure 4. Semi-unstructured Grid 
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Hanging node 


Figure 5. Single Grid Level Structure 


level 1 


level 4 


Figure 6. Multiple Grid Level Structure 





28 



face 3 


face 2 



ctof(n,l) = face 1 
ctofln^) = face 2 
ctof(n,3) = face 3 
ctof(n,4) = face 4 


Figure 7 a. Cell-to-Face Pointer 



cton(n,l) = node 1 
cton(n,2) = node 2 
cton(n,3) = node 3 
cton(n,4) = node 4 



Figure 7c. Face-to-Cell Pointer 



fton(n,l) = node 1 
ftonCn^) = node 2 



Figure 7d. Face-to-Node Pointer 


cell n 



link(n,l) = cell 1 
link(n,2) * cell 2 
link(n,3) = cell 3 
link(n,4) = cell 4 


Figure 7e. Lower Link Pointer 
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Figure 10. Inflow/Outflow Boundary 


t 



Figure 11. Characteristics at Inflow/Outflow Boundary 
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Figure 13. Computational Stencil for Fine Grid Fluxes 


Interface Boundary 



Figure 14. Computational Stencil for Coarse Grid Flux 
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Qri d L ev e l 



Figure 15. Multigrid V-cycle 
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CHAPTER IV 


GRID EMBEDDING AND ADAPTATION 

4.1 Grid and Pointer Generation 

The data structure of this algorithm requires the generation of all the pointers 
described in chapter 3. This is accomplished by first obtaining an initial structured 
grid from which all the necessary pointers are then extracted. The initial grids used 
in this study are O-type grids, as shown in figure 16. These grids are generated 
using transfinite interpolation to obtain the coordinates of each of the nodes 12 . 
Global coarser grids, required for use in the multigrid algorithm, are generated 
simply by removing every other point from the next finest grid. On each of these 
grids, the pointers are obtained in the same manner as for the finest grid. 

4.2 Adaptation 

For grid adaptation, it is first necessary to determine which regions of the 
flowfield contain high gradients. For this, various refinement parameters can be 
used which have been investigated in recent years. An extensive investigation into 
the suitability of many parameters appropriate for inviscid transonic flow has been 
conducted by Dannenlioffer 3 . A summary of the refinement parameters and their 
effectiveness in detecting flow features is given in table 4.1. As seen in the table, 
the refinement parameters |V/9| (p = denisity) and \Vq\ ( q ~ velocity magnitude) 
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Feature type 


13 

shock 

slip 

expansion 

stagnation 

shock 


wave 

line 

fan 

zone 

wake 

\p\ 

0 

0 

0 

i 

0 

VP 

2 

1 

2 

i 

0 

~ 2 
V P 

2 

2 

1 

i 

0 

IpI 

0 

0 

0 

2 - 

0 

VP 

2 

0 

2 

2 

0 

~ 2 
V P 

2 

2 

1 

2 

0 

tel 

0 

0 

0 

0 

0 ! 

vq 

2 

2 

2 

2 

0 

~ 2 
v q 

2 

2 

1 

2 

0 

bol 

1 

0 

0 

0 

2 

VPo 

2 

2 

0 

0 

2 

V 2 Po 

2 

2 

0 

0 

2 


Note: 0 => the feature is not detected 

1 => the feature is somewhat detected 

2 =>• the feature is well detected 


Table 4. 1 : Expected effectiveness of various refinement parameters for 

inviscid, transonic flows over airfoil-like bodies (from reference [3]) 
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are both capable of detecting all of the features considered except the wake region 
behind a shock. As noted in reference 3, refinement parameters based on total 
pressure can be used to successfully refine this region but tend to over-refine the 
grid downstream of the shock with the result that the grid extending from the 
shock all the way to the downstream boundary is flagged for refinement. 

The operator V used in defining the refinement parameter is the undivided 
gradient whose magnitude for an arbitrary scalar <f> is given by 

I W|= 

The partial derivatives in equation (4.1) are currently evaluated using undivided 
central differences. The refinement parameters, (3 — |Vp| and (3 — |V?|, are 
normalized to form a new parameter, p 




n 2 


+ 


d<f> 

l^J 


T 2 


(4.1) 


P = 




(3 max /3-n 


(4.2) 


which varies from 0.0 to 1.0. 

After a solution is obtained on a current grid level, each childless cell is exam- 
ined to determine whether 0 for that cell is less than a user input threshold value 
in which case the cell is flagged for refinement. Note that the value of (3* is 
somewhat arbitrary. However, selecting a threshold value too low will embed too 
many cells and generally pick up “noise” in the solution. Conversely, a value too 
high will not embed enough cells and may result in critical flow features remaining 
undetected. For the current study, it has been found that specifying (3* to be 
approximately 0.1 generally gave acceptable results. 
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After flagging cells, voids in the new grid topology can exist as shown in 
figure 17. These voids occur in areas of the grid where the refinement parameter 
j3 is close to the threshold value ,3* causing a scattering of the flagged cells over 
the parent grid. Voids in the adapted grid region can be avoided by eliminating 
undesirable cell arrangements using a mapping function along with the face-to-cell 
pointers to add new cells. This mapping function is illustrated in figure 18 in which 
five undesirable cell arrangements have been identified and a more desirable cell 
arrangement is suggested for each. This process of eliminating voids is referred to 
as grid amalgamation. 

After the grid amalgamation process, the ghost cells for the embedded inter- 
faces must be created. These cells are not constructed in the same maimer as the 
computational embedded cells since they are updated by interpolation and not 
through the time advancement scheme. Only the nodes and link pointers are gen- 
erated for these cells. The link pointers are needed for the bilinear interpolation. 
The link pointers for the parent cell do not point to the children who are embed- 
ded ghost cells since no information is passed from the children to the parents for 
these cells (such as the relative truncation error). 

After the cells are flagged for adaptation, new nodes must be generated in 
order to refine the grid in these regions. The generation of new nodes, however, 
is somewhat difficult since it is necessary to maintain both the grid smoothness 
and stretching of the parent grid. This is particularly true for cells which lie on 
the surface of the geometry since nonsmooth grid fines can result in oscillations in 
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the solution. The most straightforward approach to generating new grid points is 
to simply average the coordinates of the neighboring nodes. This technique, how- 
ever, does not resolve the curvature of the grid lines and will produce discontinous 
metrics. Another method for generating new points is to use an interpolation 
function, such as cubic splines. However, this approach presents two problems. 
The first problem arises if there is an inadequate number of defining points in the 
initial grid so that the resulting spline coefficients may not represent the true grid 
shape accurately. The second problem is caused by the pointer structure used for 
the unstructured grids which only provides local information on the connectivity 
of faces and does not provide the global connectivity required to implement an in- 
terpolation function such as cubic splines. This could be circumvented by creating 
the connectivity required but would add to the complexity of the data structure 
and would not eliminate the first problem of the initial coarse grid. 

The current approach used to construct the embedded regions is to create the 
finest level grid with a structured grid generation package and extract the nodes 
from this fine grid data file as needed. This approach, although costly, ensures 
that the new grids maintain the same stretching and smoothness. 

After the embedded regions are generated, the solution from the previous grid 
is interpolated to the new cells using the same bilinear interpolation operators 
used by the multi-grid algorithm. The solution is then advanced in time using the 
multigrid method of section 3.6. The overall adaptation process for a single grid 
level is illustrated in figure 19. 
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Figure 19. Grid Adaptation Steps 
























CHAPTER V 


RESULTS 


Results are now presented for the adaptive semi-unstructured Euler solver. 
The computer code for this algorithm used 105 storage locations for each compu- 
tational cell. The iterative scheme of the Euler solver is completely vectorizable 
and has a computational rate of 14 microseconds per cell per iteration for a sin- 
gle grid level on the Cray-2 supercomputer at NASA Langley Research Center. 
Implementation of the multigrid algorithm with five grid levels requires approxi- 
mately 30 microseconds per cell per iteration but this is dependent on the number 
of embedded ghost cells in the computational domain. All of the cases presented 
here were run using k — | for equations 3.7 and 3.8. 

Comparisons are made with the two-dimensional Euler code of reference 4 
referred to as CFL2D. As mentioned previously, the flux integration techniques 
for this code are the same as the present adaptive code. The version of CFL2D used 
for comparisons requires 139 storage locations per cell and has a computational 
rate of approximately 45 microseconds per cell per iteration. This Euler code 
uses an implicit, time advancing scheme and has better smoothing characteristics 
for the multigrid algorithm than the explicit scheme used in the adaptive code. 
Hence, on a given structured grid, CFL2D will converge in less cpu time than the 
adaptive code. 
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The grid used with CFL2D was a 257x97 O-grid with an outer boundary 
40 chords away from the airfoil. This was the same fine grid used for extracting 
points for the adaptive Euler code. The solutions presented here were considered 
converged after the rms of the average residual had been reduced by at least 
7 orders of magnitude which generally required 2000 iterations for the adaptive 
code (1 V-cycle = 1 iteration). 

5.1 NACA 0012, Mach = 0.63, a = 2.0° 

The NACA 0012 airfoil in a freestream Mach number of 0.63 and at 2.0 degrees 
angle-of-attack is presented here. This is an isentropic case and should produce 
a drag coefficient of zero. Figure 20 shows the initial grid for the adaptive Euler 
code which is a 33x13 structured O-grid. The surface pressure coefficients are 
shown in figure 21 and the Mach number contours are shown in figures 22 and 23. 
The coarseness of the grid produces large entropy values near the leading edge of 
the airfoil and degrades the solution downstream which is evident from the poor 
resolution of the trailing edge stagnation point. 

The first adaptive embedding was performed using the undivided gradient 
of the density magnitude | V p | with the threshold value p* set to 0.1. The 
solution was converged on each level before the next level of embedded grids was 
created. Figures 24 through 29 show each successive level of refinement with the 
corresponding surface pressures. As shown, the embedding significantly improves 
the overall solution quality without a large increase in the number of cells. Figures 
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36 and 37 show the Mach contours on the final adapted grid. These contours show 
the smooth transition of the solution at the embedded interfaces. 

The second adaptive embedding case was performed using the undivided gra- 
dient of the velocity magnitude \ Vq \ with f3* again set to 0.1. The computations 
were carried out in the same manner as the previous case with the resulting grids 
and surface pressures shown in figures 30 through 35. As shown, this refinement 
parameter embedded more points than the density gradient for the same value of 
/?*. Table 5.1 shows a comparison of the lift and drag values obtained for the two 
cases along with those obtained with CFL2D. 

5.2 NACA 0012, Mach = 0.8, a = 1.25° 

The NACA 0012 airfoil in a frees tream Macli number of 0.80 and at 1.25 
degrees angle of attack is presented here. Using | V p | or | Vq | produced nearly 
identical grids for this case and so only the | Vg| case is shown. 

These freestream conditions produce a strong shock on the upper surface of 
the airfoil with a relatively weak shock on the lower surface. The Mach contours 
and pressure distribution for the initial grid are shown in figures 38 and 39. The 
final adapted grid, pressure distribution, and Mach contours are shown in figures 
40 through 42. As shown, the final grid produces a sharp shock on both Hie 
upper and lower surfaces with no oscillations. Table 5.2 shows the comparison 
with CFL2D. The final c* has about a 4 ^ % difference from that of CFL2D. 
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Figure 29. 4 Level (Vp) Pressure Dis 
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Adaptation 

Criteria 

Number of 
Cells 

q 

c d 

None 

(CFL2D) 

32,298 

.328 

0.270 x 10' 3 

Vq 

11,132 

.320 

0.249 x 10‘ 3 

Vp 

7,648 

.323 

0.769 x 10' 3 


Table 5.1. Results Comparison for M«, = 0.63, a = 2.0° 
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Figure 42. 4 Level (Vq) Mach Contours, Moo = 0.80, a - 1.25° 


Adaptation 

Criteria 

Number of 
Cells 


c d 

None 

(CFL2D) 

32,298 

.352 

0.211 x 10' 1 

Vq 

11,132 

.337 

0.251 x 10' 1 


Table 5.2. Results Comparison for Moo = 0.80, a = 1.25° 
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CHAPTER VI 


CONCLUDING REMARKS 

A semi-unstructured algorithm to solve the two-dimensional Euler equations 
with adaptive grid embedding has been presented, A multigrid algorithm has been 
implemented to provide the required grid level communication and to accelerate 
the convergence. 

Results have been shown for a subcritical airfoil and comparisons have been 
made with a structured algorithm using the same flux integration techniques. 
The final lift coefficients are within 2|% of those obtained with CFL2D. Results 
have also been shown for a transonic airfoil where the increased resolution of the 
shocks is very substantial without the large increase in grid points. The final lift 
coefficients for this case are around 4|% of CFL2D. 

One problem with the present method is the poor error-damping character- 
istics of the two-stage time advancing scheme. Further research is needed in this 
area to produce an acceptable scheme. The memory requirements on a per cell 
basis for the present method are approximately the same as that for an implicit 
scheme. Since the present method requires less grid points to resolve the flow, the 
overall memory requirements are less. The advantages of using this scheme would 
probably be greater in three-dimensions. 
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APPENDIX A 


Interpolation Operators 


Interpolation can be performed in a two-dimensional field, shown in figure 43, 
by assuming a bilinear function of the form: 


/(*, y) = a + bz + cy + dxy 


(A.l) 


The constants a, 6, c, and d, can be evaluated given four boundary values: 


= a + bxi + cyi + dxij/i (.4.2a) 

f{x2,y2) = a + b*2 + W2 + dz 2 y2 {A.2b) 

f(x 3 , y 3 ) = a + bx 3 + cy 3 + dx 3 y 3 {A. 2c) 

/(z-hJAi) = a + bx 4 + cy 4 + dx A y A (A. 2d) 


This can be written in matrix form as: 

( i *1 yi *iyi\ /«\ //(*i»yi)\ 

1 *2 V2 *2V2 | I b I I /(*2,y2) I i A 0\ 

l *3 y 3 * 3 y 3 I I c I - I /(* 3 ,y 3 ) I 1 ' 

1 *4 1/4 X 4 y 4 / Vd/ V/(^4,y4)/ 

Tliis matrix can be simplified if (aci,yi) is defined to be (0,0). Equation (A. 2a) is 
reduced to: 


f{*uyi) = a M- 4 ) 

This reduces (A.3) to a 3x3 matrix. 

( a-2 y 3 *2V2 \ / b\ / A f 2 \ 

3\ 3 y 3 * 3 y 3 ) I c ) = I A/3 J (d.. 5 ) 

*4 y4 3 - 4 y 4 J \dj \ A/4 J 
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where: 


A/ n = f{x n ,y n ) - /(*i,y l) 


The system of equations given in (A. 5) can now easily be solved by the use of 
Kramer’s rule. The determinant of the coefficient matrix in (^4.5) is given as: 


A = 


*2 

y2 

*2y2 

*3 

ys 

a-aya 

^4 

y4 

a-4y4 


*2*41/3^4 -y2) + *2*3y4(y2 -y3) + *3*4y2(y3 -y4) (-4-6) 


The coefficients can then be solved for: 

A/ 2 y 2 x 2 y 2 
Af 3 y3 x 3 y 3 
A/ 4 y 4 x 4 y 4 


6 = 


y3y*{x4 - ^3)A/ 2 + y 2 y4(»2 - s4)A/ 3 + y 2 y 3 (;c3 - s2)A/ 4 

A 


(A.7a) 


c — 


^2 

a/ 2 

^2y2 

3*3 

a/ 3 

*3y3 

*4 

a/ 4 



X 3 X 4 (y 3 - y 4 )A/ 2 + X 2 X 4 (y 4 - V 2 )A/ 3 + X 2 X 3 (y2 - y 3 )A/ 4 


(A.76) 


*2 

y 2 

a/ 2 

3-3 

ys 

a/ 3 

*4 

y4 

a/ 4 


(a: 3 y 4 - x 4 y 3 )Af 2 + (z 4 y 2 - x 2 y 4 )Af 3 + (z 2 y 3 - x 3 y 2 )Af 4 


(A.7c) 


Assuming a uniform Cartesian grid of unit width and lieighth, the nodes would 

(*i,yi) = (0,0) 

(x 2 ,y 2 ) = (1,0) 

(*3,ys) = (1,1) 


(-4.8) 


(a-4,y4) = (0,1) 
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Substituting (>1.8) into (>1.7) yields the following solution: 


a = /(*i,yi) 

b = A/ 2 

(>L.9) 

c = Af 4 

d=- A/ 2 + A/ 3 - A/ 4 
Substituting (>1.9) into (>1.1) yields: 

/(*%y) =/(*i,yi) + (/(* 2 ,y 2 ) - /(*i,yi))* + (/(*4,y4) - /(*i,yi))y 

+ (-(/(*2 ? y2) - /(*i»vi)) (>1.10) 

+ (/(*3,ys) - f(xi,yi)) - {f{x 4 ,y 4 ) - f{x u yi)))xy 
Simplifying this gives: 

f{x,y) = (1 - x - y + xy)f(x 1 ,y 1 ) + (x - xy)f(x 2 ,y 2 ) 

(.4.11) 

+ (*y)/(*3,y3) + (y - *y)/(* 4,y4) 

For a Cartesian grid as shown in figure 44, (x,y) = (0.25,0.25), so that 
/(0.25, 0.25) = (0.5625)/(ar 1 ,y 1 ) + (0.1875)/(* 2 ,y 2 ) 

(> 1 . 12 ) 

+ (0.0625)/(* 3 ,y 3 ) + (0.1875)/(a- 4 ,t/ 4 ) 

Although the grids are generally not Cartesian, this is a close approximation to 
the actual weightings for grids that are not highly stretched. 
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Figure 43. Typical Quadrilateral 



• - Coarse grid cell-centered data 


O - Fine grid cell-centered data 


Figure 44. Quadrilateral in a Uniform Cartesian Grid 
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